# misc.fun_trans_elem -- 
#
# Copyright (c) 2001, Willem Jan Zaadnoordijk

from math        import exp,pi,sqrt

from geometry    import XY
from misc.funE_B import E1,erfc
from misc.funLB  import ieer


def wpp1(xyw,tt,rp):
	'''function for unit Theis well potential'''
	return (-0.25/pi*E1( .25*(xyw.x**2+xyw.y**2) / (tt/rp) ))


def slp1(halflength, czl, rtime, rp) :
	"potential @czl for constant step-linesink rtime after start"
#  potential at CZ, RTIME for line-sink from (-halflength,0 to halflength,0)CZ1 to CZ2 with
#  constant strength both in time and in space starting at RTIME=0
#  with transient constant RP
#  (=porosity/(hydaulic conductivity*average head
	rsqrpi=1.77245385090551
	rtolrl=.0001
	rrfar =2.
#  maybe test for very small (dimensionless) time in stead of zero
	if (rtime <= 0.) :
		return 0
#  the element is translated and rotated such that the center is at
#  the origin and the element along the x-axis with endpoint 1 at x1=
#  -.5*(length of element), endpoint 2 at x2=.5*(length of element)
	rl  = halflength*2
	rdum= abs(czl)/rl
	if (rdum > rrfar) :
		# far-field function: unit well times length
		rfslp1 = wpp1(XY(czl.real,czl.imag), rtime, rp) * rl
	else :
		rx1=-halflength
		rx2=-rx1
		rtol=rtolrl*2.*rx2
		rx=czl.real
		ry=czl.imag
		if (abs(czl-rx1) < rtol) :
			rx=rx1-rtol
		elif (abs(czl-rx2) < rtol) :
			rx=rx2+rtol
		if (abs(ry) < rtol) :
			if (ry < 0.) :
				ry=-rtol
			else :
				ry=rtol
		rpo4t =.25*rp/rtime
		rsqpo4=sqrt(rpo4t)
		rlim=rsqpo4*abs(ry)
		r1=-sqrt(rtime/rp)/rsqrpi*.5*exp(-rpo4t*ry*ry)*   \
			( erfc(rsqpo4*(rx-rx2))-erfc(rsqpo4*(rx-rx1)) )
		r2=.25/pi*   \
			( (rx-rx2)*E1(rpo4t*((rx-rx2)*(rx-rx2)+ry*ry))   \
			 -(rx-rx1)*E1(rpo4t*((rx-rx1)*(rx-rx1)+ry*ry)) )
		r3=ry*.5/rsqrpi*ieer( rlim, (rx-rx2)/ry, (rx-rx1)/ry )
		rfslp1=r1+r2+r3
	return rfslp1


def slq1(halflength, czl, rtime, rp) :
#  Qx and Qy at CZ, RTIME for line-sink from CZ1 to CZ2 with
#  constant strength both in time and in space starting at RTIME=0
#  with transient constant RP
#  (=porosity/(hydaulic conductivity*average head
	rsqrpi=1.77245385090551
	rtolrl=.0001
#  maybe test for very small (dimensionless) time in stead of zero
	if (rtime <= 0.) :
		cfslq1=XY(0.,0.)
		return cfslq1
#  the element is translated and rotated such that the center is at
#  the origin and the element along the x-axis with endpoint 1 at x1=
#  -.5*(length of element), endpoint 2 at x2=.5*(length of element)
	rx1=-halflength
	rx2= halflength
	rtol=rtolrl*2.*rx2
	rx=czl.real
	ry=czl.imag
	if (abs(czl-rx1) < rtol) :
		rx=rx1-rtol
	elif (abs(czl-rx2) < rtol) :
		rx=rx2+rtol
	if (abs(ry) < rtol) :
		if (ry < 0.) :
			ry=-rtol
		else :
			ry=rtol
	rpo4t =.25*rp/rtime
	rlim=sqrt(.25*rp/rtime)*ry
	rqxi=-.25/pi*( E1(rpo4t*((rx-rx2)*(rx-rx2)+ry*ry))   \
	              -E1(rpo4t*((rx-rx1)*(rx-rx1)+ry*ry)) )
	rqeta=-.5/rsqrpi*ieer( rlim, (rx-rx2)/ry, (rx-rx1)/ry )
	cfslq1=XY(rqxi,rqeta)
	return cfslq1
